#ifndef _XMESH_EXACT_H
#define _XMESH_EXACT_H
#include <cstdio>
#include <cstdlib>
#include <typeinfo>

#define REAL double 

#define Fast_Two_Sum_Tail(a, b, x, y)			\
	bvirt = x - a;								\
	y = b - bvirt

#define Fast_Two_Sum(a, b, x, y)				\
	x = (REAL) (a + b);							\
	Fast_Two_Sum_Tail(a, b, x, y)

#define Fast_Two_Diff_Tail(a, b, x, y)			\
	bvirt = a - x;								\
	y = bvirt - b

#define Fast_Two_Diff(a, b, x, y)				\
	x = (REAL) (a - b);							\
	Fast_Two_Diff_Tail(a, b, x, y)

#define Two_Sum_Tail(a, b, x, y)				\
	bvirt = (REAL) (x - a);						\
	avirt = x - bvirt;							\
	bround = b - bvirt;							\
	around = a - avirt;							\
	y = around + bround

#define Two_Sum(a, b, x, y)						\
	x = (REAL) (a + b);							\
	Two_Sum_Tail(a, b, x, y)

#define Two_Diff_Tail(a, b, x, y)				\
	bvirt = (REAL) (a - x);						\
	avirt = x + bvirt;							\
	bround = bvirt - b;							\
	around = a - avirt;							\
	y = around + bround

#define Two_Diff(a, b, x, y)					\
	x = (REAL) (a - b);							\
	Two_Diff_Tail(a, b, x, y)

#define Split(a, ahi, alo)						\
	c = (REAL) (splitter * a);					\
	abig = (REAL) (c - a);						\
	ahi = c - abig;								\
	alo = a - ahi

#define Two_Product_Tail(a, b, x, y)			\
	Split(a, ahi, alo);							\
	Split(b, bhi, blo);							\
	err1 = x - (ahi * bhi);						\
	err2 = err1 - (alo * bhi);					\
	err3 = err2 - (ahi * blo);					\
	y = (alo * blo) - err3

#define Two_Product(a, b, x, y)					\
	x = (REAL) (a * b);							\
	Two_Product_Tail(a, b, x, y)

/* Two_Product_Presplit() is Two_Product() where one of the inputs has       */
/*   already been split.  Avoids redundant splitting.                        */

#define Two_Product_Presplit(a, b, bhi, blo, x, y)	\
	x = (REAL) (a * b);								\
	Split(a, ahi, alo);								\
	err1 = x - (ahi * bhi);							\
	err2 = err1 - (alo * bhi);						\
	err3 = err2 - (ahi * blo);						\
	y = (alo * blo) - err3

/* Two_Product_2Presplit() is Two_Product() where both of the inputs have    */
/*   already been split.  Avoids redundant splitting.                        */

#define Two_Product_2Presplit(a, ahi, alo, b, bhi, blo, x, y)	\
	x = (REAL) (a * b);											\
	err1 = x - (ahi * bhi);										\
	err2 = err1 - (alo * bhi);									\
	err3 = err2 - (ahi * blo);									\
	y = (alo * blo) - err3

/* Square() can be done more quickly than Two_Product().                     */

#define Square_Tail(a, x, y)					\
	Split(a, ahi, alo);							\
	err1 = x - (ahi * ahi);						\
	err3 = err1 - ((ahi + ahi) * alo);			\
	y = (alo * alo) - err3

#define Square(a, x, y)							\
	x = (REAL) (a * a);							\
	Square_Tail(a, x, y)

/* Macros for summing expansions of various fixed lengths.  These are all    */
/*   unrolled versions of Expansion_Sum().                                   */

#define Two_One_Sum(a1, a0, b, x2, x1, x0)		\
	Two_Sum(a0, b , _i, x0);					\
	Two_Sum(a1, _i, x2, x1)

#define Two_One_Diff(a1, a0, b, x2, x1, x0)		\
	Two_Diff(a0, b , _i, x0);					\
	Two_Sum( a1, _i, x2, x1)

#define Two_Two_Sum(a1, a0, b1, b0, x3, x2, x1, x0) \
	Two_One_Sum(a1, a0, b0, _j, _0, x0);			\
	Two_One_Sum(_j, _0, b1, x3, x2, x1)

#define Two_Two_Diff(a1, a0, b1, b0, x3, x2, x1, x0)	\
	Two_One_Diff(a1, a0, b0, _j, _0, x0);				\
	Two_One_Diff(_j, _0, b1, x3, x2, x1)

#define Four_One_Sum(a3, a2, a1, a0, b, x4, x3, x2, x1, x0) \
	Two_One_Sum(a1, a0, b , _j, x1, x0);					\
	Two_One_Sum(a3, a2, _j, x4, x3, x2)

#define Four_Two_Sum(a3, a2, a1, a0, b1, b0, x5, x4, x3, x2, x1, x0)	\
	Four_One_Sum(a3, a2, a1, a0, b0, _k, _2, _1, _0, x0);				\
	Four_One_Sum(_k, _2, _1, _0, b1, x5, x4, x3, x2, x1)

#define Four_Four_Sum(a3, a2, a1, a0, b4, b3, b1, b0, x7, x6, x5, x4, x3, x2, \
                      x1, x0)											\
	Four_Two_Sum(a3, a2, a1, a0, b1, b0, _l, _2, _1, _0, x1, x0);		\
	Four_Two_Sum(_l, _2, _1, _0, b4, b3, x7, x6, x5, x4, x3, x2)

#define Eight_One_Sum(a7, a6, a5, a4, a3, a2, a1, a0, b, x8, x7, x6, x5, x4, \
                      x3, x2, x1, x0)									\
	Four_One_Sum(a3, a2, a1, a0, b , _j, x3, x2, x1, x0);				\
	Four_One_Sum(a7, a6, a5, a4, _j, x8, x7, x6, x5, x4)

#define Eight_Two_Sum(a7, a6, a5, a4, a3, a2, a1, a0, b1, b0, x9, x8, x7, \
                      x6, x5, x4, x3, x2, x1, x0)						\
	Eight_One_Sum(a7, a6, a5, a4, a3, a2, a1, a0, b0, _k, _6, _5, _4, _3, _2, \
				  _1, _0, x0);											\
	Eight_One_Sum(_k, _6, _5, _4, _3, _2, _1, _0, b1, x9, x8, x7, x6, x5, x4, \
				  x3, x2, x1)

#define Eight_Four_Sum(a7, a6, a5, a4, a3, a2, a1, a0, b4, b3, b1, b0, x11, \
                       x10, x9, x8, x7, x6, x5, x4, x3, x2, x1, x0)		\
	Eight_Two_Sum(a7, a6, a5, a4, a3, a2, a1, a0, b1, b0, _l, _6, _5, _4, _3, \
				  _2, _1, _0, x1, x0);									\
	Eight_Two_Sum(_l, _6, _5, _4, _3, _2, _1, _0, b4, b3, x11, x10, x9, x8, \
				  x7, x6, x5, x4, x3, x2)

/* Macros for multiplying expansions of various fixed lengths.               */

#define Two_One_Product(a1, a0, b, x3, x2, x1, x0)	\
	Split(b, bhi, blo);								\
	Two_Product_Presplit(a0, b, bhi, blo, _i, x0);	\
	Two_Product_Presplit(a1, b, bhi, blo, _j, _0);	\
	Two_Sum(_i, _0, _k, x1);						\
	Fast_Two_Sum(_j, _k, x3, x2)

#define Four_One_Product(a3, a2, a1, a0, b, x7, x6, x5, x4, x3, x2, x1, x0) \
	Split(b, bhi, blo);													\
	Two_Product_Presplit(a0, b, bhi, blo, _i, x0);						\
	Two_Product_Presplit(a1, b, bhi, blo, _j, _0);						\
	Two_Sum(_i, _0, _k, x1);											\
	Fast_Two_Sum(_j, _k, _i, x2);										\
	Two_Product_Presplit(a2, b, bhi, blo, _j, _0);						\
	Two_Sum(_i, _0, _k, x3);											\
	Fast_Two_Sum(_j, _k, _i, x4);										\
	Two_Product_Presplit(a3, b, bhi, blo, _j, _0);						\
	Two_Sum(_i, _0, _k, x5);											\
	Fast_Two_Sum(_j, _k, x7, x6)

#define Two_Two_Product(a1, a0, b1, b0, x7, x6, x5, x4, x3, x2, x1, x0) \
	Split(a0, a0hi, a0lo);												\
	Split(b0, bhi, blo);												\
	Two_Product_2Presplit(a0, a0hi, a0lo, b0, bhi, blo, _i, x0);		\
	Split(a1, a1hi, a1lo);												\
	Two_Product_2Presplit(a1, a1hi, a1lo, b0, bhi, blo, _j, _0);		\
	Two_Sum(_i, _0, _k, _1);											\
	Fast_Two_Sum(_j, _k, _l, _2);										\
	Split(b1, bhi, blo);												\
	Two_Product_2Presplit(a0, a0hi, a0lo, b1, bhi, blo, _i, _0);		\
	Two_Sum(_1, _0, _k, x1);											\
	Two_Sum(_2, _k, _j, _1);											\
	Two_Sum(_l, _j, _m, _2);											\
	Two_Product_2Presplit(a1, a1hi, a1lo, b1, bhi, blo, _j, _0);		\
	Two_Sum(_i, _0, _n, _0);											\
	Two_Sum(_1, _0, _i, x2);											\
	Two_Sum(_2, _i, _k, _1);											\
	Two_Sum(_m, _k, _l, _2);											\
	Two_Sum(_j, _n, _k, _0);											\
	Two_Sum(_1, _0, _j, x3);											\
	Two_Sum(_2, _j, _i, _1);											\
	Two_Sum(_l, _i, _m, _2);											\
	Two_Sum(_1, _k, _i, x4);											\
	Two_Sum(_2, _i, _k, x5);											\
	Two_Sum(_m, _k, x7, x6)

/* An expansion of length two can be squared more quickly than finding the   */
/*   product of two different expansions of length two, and the result is    */
/*   guaranteed to have no more than six (rather than eight) components.     */

#define Two_Square(a1, a0, x5, x4, x3, x2, x1, x0)	\
	Square(a0, _j, x0);								\
	_0 = a0 + a0;									\
	Two_Product(a1, _0, _k, _1);					\
	Two_One_Sum(_k, _1, _j, _l, _2, x1);			\
	Square(a1, _j, _1);								\
	Two_Two_Sum(_j, _1, _l, _2, x5, x4, x3, x2)


REAL splitter ;
REAL epsilon ;

static void printdouble (double number)
{
	unsigned long long no;
	unsigned long long sign, expo;
	int exponent;
	int i, bottomi;

	no = *(unsigned long long *) &number;
	sign = no & 0x8000000000000000ll;
	expo = (no >> 52) & 0x7ffll;
	exponent = (int) expo;
	exponent = exponent - 1023;
	if (sign) {
		printf("-");
	} else {
		printf(" ");
	}
	if (exponent == -1023) {
		printf(
			"0.0000000000000000000000000000000000000000000000000000_     (   )");
	} else {
		printf("1.");
		bottomi = -1;
		for (i = 0; i < 52; i++) {
			if (no & 0x0008000000000000ll) {
				printf("1");
				bottomi = i;
			} else {
				printf("0");
			}
			no <<= 1;
		}
		printf("_%d  (%d)", exponent, exponent - 1 - bottomi);
	}
}

void expansion_print(int elen, REAL * e)
{
	int i;

	for (i = elen - 1; i >= 0; i--) {
		printdouble(e[i]);
		if (i > 0) {
			printf(" +\n");
		} else {
			printf("\n");
		}
	}
}
void exactinit ()
{
	REAL half;
	REAL check, lastcheck;
	int every_other;

	every_other = 1;
	half = 0.5;
	epsilon = 1.0;
	splitter = 1.0;
	check = 1.0;
	/* Repeatedly divide `epsilon' by two until it is too small to add to    */
	/*   one without causing roundoff.  (Also check if the sum is equal to   */
	/*   the previous sum, for machines that round up instead of using exact */
	/*   rounding.  Not that this library will work on such machines anyway. */
	do {
		lastcheck = check;
		epsilon *= half;
		if (every_other) {
			splitter *= 2.0;
		}
		every_other = !every_other;
		check = 1.0 + epsilon;
	} while ((check != 1.0) && (check != lastcheck));
	splitter += 1.0;
}

static void printfloat (float number )
{
	unsigned no;
	unsigned sign, expo;
	int exponent;
	int i, bottomi;

	no = *(unsigned *) &number;
	sign = no & 0x80000000;
	expo = (no >> 23) & 0xff;
	exponent = (int) expo;
	exponent = exponent - 127;
	if (sign) {
		printf("-");
	} else {
		printf(" ");
	}
	if (exponent == -127) {
		printf("0.00000000000000000000000_     (   )");
	} else {
		printf("1.");
		bottomi = -1;
		for (i = 0; i < 23; i++) {
			if (no & 0x00400000) {
				printf("1");
				bottomi = i;
			} else {
				printf("0");
			}
			no <<= 1;
		}
		printf("_%3d  (%3d)", exponent, exponent - 1 - bottomi);
	}
}


template <typename ValueType>
class Real
{
public:
	typedef ValueType VT ;
public:
	Real (): x(ValueType(0)), y(ValueType(0)) {}
	Real (ValueType num) 
		: x(num), y (ValueType(0))
	{
	}

Real(const Real & _d) :  x(_d.x), y (_d.y)
	{
		
	}

	
	friend Real operator + (const Real & lhs, const Real & rhs)
	{
		
	}

	friend std::ostream & operator << (std::ostream & ostr, const Real & _n )
	{
		
		print (_n.x) ; printf ("\n") ;
		print (_n.y) ; printf ("\n") ;
		return ostr ;
	}
public :
	Real two_sum (ValueType a, ValueType b)
	{
		Real r ;
		r.x = a + b ;
		double bv = r.x - a;
		double av = r.x - bv ;
		double br = b - bv ;
		double ar = a - av ;
		r.y = ar + br ;
		return r;
	}

	static void print (ValueType num)
	{
		if (typeid (ValueType) == typeid (double))
		{
			printdouble (num );
		}
		else if (typeid (ValueType) == typeid (float))
		{
			printfloat (num );
		}
		else 
		{
			throw ("The internal type is not supported. ");
		}
	}
private :


public :
	
private :
	ValueType x, y ;
};
#endif
